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Abstract 

Recurrence plot is a quite easy tool to be used in time series analysis, in particular for measuring 
unstable periodic orbits embedded in a chaotic dynamical system. Recurrence quantified analysis 
(RQA) is an advance tool that allows the study of intrinsic complexity of dynamical system with a 
set of few parameters. We use RQA for measuring chaotic transitions of NADH chemical reaction 
and determine numerically its characteristic parameters such as: Correlation integral, information 
entropy, Maximal Lyapunov's exponent, etc. For this work we have developed command sets with 
performance better than TISEAN package 
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I. INTRODUCTION 



A. Recurrence quantified analysis 

Recurrence Plot (RP) was initially introduced by Eckman et al. (1987) [1[ as a tool for 
analyzing experimental time series data, especially useful for finding hidden correlations in 
highly complicated data and determine the stationarity of the time series. This method 
allows the identification of system properties that cannot be observed by the linear and 
nonlinear usual approaches It is worth to mention the simplicity of the algorithms during 
numerical calculations too. A RP is an injective application of a single reconstructed tra- 
jectory to the boolean matrix space, each pair y iy yj coming from the time series is related 
with a pair (i, j), called recurrence points. Let us consider N values of a time series given 
by 

{x , • • • , Xn-i}, 

with N large enough in order to evaluate the embedding dimension by using the false 
nearest neighborhood (d > 2) and the time delay (r > 1) by looking at the relative mini- 
mum in the mutual information || Following Takens' embedding theorem [[| , the dynamics 
can be appropriately represented by the phase space trajectory reconstructed by using the 
time delay vectors 

Hi — (Xi, Xi+ T , • • • , Xj-f(jV-l) r)- 

and the recurrence matrix is: 

R (i,j) = ®($h -\\Vi- UjWoo) - @(<$« -\\Vi- UjWoo) (1) 

where G is the Heaviside function and the matrix is symmetric. That means a RP is 
built by comparing all delayed vectors with each other. A dark dot is plotted, [Ruj) = 1) 
with integer coordinates when 5i < \\y% — UjWoo < $h, otherwise a white dot is plotted 
(Ruj) = 0). The interval [5i, 5h] it is known as threshold corridor. The choice of this interval 
is critical, too large produces a saturation of the RP including irrelevant points, and too 
narrow losses information. Since up to now in the literature there are not a satisfactory 
solution, an educated guess should be appropriate. In this work we used [cr/10 5 , a/10 2 ], 
where a is the standard deviation. 
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Webber et al. in order to characterize and analyze recurrent plots introduced a set of 
quantifiers, which are collectively called recurrence quantified analysis (RQA). The first of 
this quantifiers is the % recurrence (%REC), defined as: 

N 

°/ REC = 100 (2) 
where N t = dim(i?) (every possible points) and N r is number of recursive points given 

by: 

N r = 2#{(i,j)/R iid) > Oandt < j} 

The slope of the linear region in the S-shaped %REC vs. corridor width is the correlation 
dimension. The second RQA quantifier is called % determinism (%DET); and it is defined 
as: 

°/„DET = 100 j-f- (3) 
where Ni is called the number of periodic dots given by: 

Ni = 2#{(i,j)/(i,j) e d c (k,b), % < j, Vc, k, b > 0} 
and a periodic line with length b, origin k and zone c > is defined as: 

k+b 

d c (k, 6) = {(i, i + c)/ J] R(i,i+c) > 0} 

i=k 

The %DET is related with the organization of the RP. The third RQA quantifier, called 
entropy (S), is closely related to %DET. 

S = -X>log 2 (P fe ) (4) 

6=1 

where H is the length of the maximum periodic line, P n ^ is the relative frequency 
of the periodic lines with length b > 0. The label entropy is just that, a label, not to be 
confused with Shannon's entropy since there is not a one to one correspondence between 
this quantifier and the Shannon's entropy. This quantity should be labeled more properly 
as first rate cumulant since is related with the relative frequency fluctuations. Moreover, 
for periodic orbits they are mapped onto diagonals with different lengths and uniformly 
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distributed, giving values of S > 0, with a maximum value. Webber assumes that S is 
related with Shannon's entropy if and only if the system is chaotic and the embedding 
dimension large enough. The fourth quantifier is the longest periodic line found during the 
computation of %DET given by L max . Eckman et al. claim that line lengths on RP are 
directly related to inverse of the largest positive Lyapunov exponent. Short lines values are 
therefore indicative of chaotic or stochastic behavior. 



B. Dynamical system, a description 

The peroxidas-oxidase (PO) reaction is an important bridge between the chemical ex- 
citable oscillators Beluzov Zabotinsky reaction and biological oscillators such as intracellular 
Ca 2+ oscillators 0. It is now clear that PO reactions shows a wide spectrum of interest- 
ing behaviors including simple oscillations, bistability, quasiperiodicity and chaos [|7| . These 
behavior have all been observed in vitro under well controlled laboratory conditions. This re- 
action appears in plants as part of the process of lignifications || with nicotinamide adenine 
dinucleotide (NADH) as electron donor. Its estequeometric formulae is: 

2 NADH + 2 + 2 H -> 2 NAD+ + 2 H 2 

In 1983, a model of PO reaction, now commonly referred to the Olsen's Model ||, was 
proposed. Simulations with the Olsen model quantitatively reproduce both the simple and 
chaotic oscillations of PO reactions. Studies with this model showed that increasing the 
parameter k 3 caused the system to undergo a transition from simple oscillations to chaos 
via a cascade of periodic doubling bifurcations. The Olsen's model involves four variable, 
molecular oxygen A, NADH B and two intermediate species. One of the intermediates is 
very likely NAD + X, while other is oxyferrous peroxidase Y, also know as compound III. 
The complete mechanism corresponds system of four differential equations are given below: 



dA 

= k 7 - k_ 7 A - k 3 ABY (5) 

^ = k 8 - k x BX - k 3 ABY (6) 
dX 

— = k x BX - 2k 2 X 2 + 3k 3 ABY - k 4 X + k 6 (7) 
dY 

— = 2k 2 X 2 - k 5 Y - k 3 ABY (8) 
dt 



We use ks as control parameter which induces a transition to chaos type I 
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II. EXPERIMENTS AND RESULTS 
A. Experimental Setup 

If we wisht to use the RQA method to a time series long enough for reliability, we face 
with the usual difficulties, memory capacity and calculus complexity Most of the codes 
were developed in C under Linux operating system. This language does not implement the 
Boolean data type using integer type instead, so the required memory amount to store the 
recurrence matrix coming from a typical time series of 2 16 data points is 128 Gbits and, the 
algorithm complexity is quadratic, which implies huge time consuming. In order to solve 
these problems we propose to reduce the dimension of the matrix R, but then the mapping 
may not be one to one and a whole region of the reconstructed space would have the same 
points as an image. Therefore R as defined in ([!]) is modified as: R(k,i) > if and only if at 
least a pair of indexes exists 

(i,j)£{(hj)/k = [i-],l = \j^)} 
m m 

for k y I given, with N = dim(i?) and m the amount data points, such as 
Q(S h - Wyi - VjWoo) - Q(5i -\\yi- VjWoo) > 

and zero otherwise. 

Due to hardware limitation we use recurrence matrices 10240 x 10240 values, although 
this fact was discussed before, and means a loss in the matrix resolution, the obtained results 
are quite satisfactory taking into account an adjusting of the corridor width and studying 
the parameters asymptotic behavior as a function of the number of the data. 

Another problem concers the RQA high sensibility with the number of cycles during a time 
interval. As is depicted in figure-|l] %DET has an asymptotic value when the number of cycles 
is greater than 30. If we are dealing with a sub-harmonic cascade or in the intermittency 
region an abrupt change in the number of cycles is observed, giving us false transitions 
associate to the change in the cycles number. But the number of cycles can not be so big 
because produces a graphics saturation then %DET falls on minor values. A reason is an 
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increment in vectors amount mapped to same (k, I) matrix component. So when this amount 
goes too far a heuristic limit; determinism is lost. Taking this into account the parametric 
measurements associate to the RQA were done in such a way to use the half cycles on the 
plateau. 




FIG. 1: %DET as a function of cycles number for PO chemical reaction. For k% = 0.028 a limit 
cycle periodic behavior is observed 



B. Application to PO reaction 

Since RP is not quite sensible to the election of the embedding dimension, d, we use 
d = 9 which this value allow us to eliminate the non diagonal points in the RP plots without 



affecting its general structure ||11|| . On the other hand, no difference is observed neither 
S(k 3 ), or any RQA parameters within a range 3 > d > 15, so we use the value which 
recommended by Takens' theorem |4]]. As can be seen in figures-^], |3], and |]it was no found 
any direct proportion between S(ks), and the maximum Lyapunov exponent, X(ks), by using 
the Aurell's algorithm et al. in TISEAN package [113] calculated with a Thieler's window 
| 13| with 500 time units. Different algorithms associate with the evaluation of X(k 3 ) were 
used with unsatisfactory results due to the time consuming for getting a reliable slope in 
the graphical methods, so we developed our own software that allow us to apply the RQA 
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methods for the calculation of the relevant parameters which characterize in any dynamical 
system independently of the size of the data. It is important to note as depicted in figure-^ 
the chaos transition characterization since S(kz) increase its value towards a plateau for the 
&3 values which correspond to a sub-harmonic cascade, as well as a peak in the neighborhood 
of the first bifurcation value. 
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FIG. 2: S(ks) for 40 cycles, a plateau can be seen as in the chaotic as in the periodic regions. 



7 




FIG. 3: Maximum Lyapunov exponent, \(ks), found by using Aurell's et al. algorithm for 40 
cycles. 




FIG. 4: 7 (fc 3 ) = log 2 (imax) for 40 cycles. 

%DET(/c 3 ) is associated with the existence of stable and unstable periodic orbits. In 
general the probability density related with periodic orbits is uniform, but for a dense 



8 



family of unstable it orbits has an exponential distribution since it is more probable to find 
our systems orbits with shorter period. For us, %DET>99% always, see figure-|5], and the 
behavior in the chaotic region is the same as S(ks), that means the increasing number of 
periodic orbits due to the sub-harmonic cascade from fc 3 = 0.031 up to reach a maximum 
in the chaotic region densely populated by unstable orbits. For k 3 > 0.0355 the number 
of orbits is considerably reduced after an intermittent behavior. This is intrinsic to our 
dynamical system since it is not observed in the Lorentz system, used as a toy model |Tl|]in 
this trail. 
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FIG. 5: %DET as a function of k^,it is not observed relevant changes in the curves shape due to 
the abundance of unstable periodic orbits. 

On the other hand, by analyzing %KEC(k 3 ) we may appreciate the scarce population of 
recurrent points in the RP. %REC holds a constant value and is reduced at the beginning of 
the cascade with a minimum in the chaotic region. This behavior is different of that observed 
with %DET and S augmenting the value up to get a lower plateau. We assume this can be 
due to the changes in the shape of the attractor, since the same results is obtained for the 



Lorentz system |L1 . 
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FIG. 6: %REC as a function of ks,it can be seen how scarcely populated is RP. 

By a simple inspection in the RP the increasing number of periodic orbits is clear, because 
more parallel lines to the diagonal appear, but the changes in the attractor shape are less 
impressive and only can be appreciate as alabeated curves when the embedding dimension 
is near d = 3. 

III. CONCLUSIONS 

In spite of Webber's assertion we do not observed any difference in the S(k 3 ) graph 
neither for high embedding dimension, d > 9, nor the lower ones, 2 < d < 10. S(k 3 ) shows a 
direct proportionality with the maximum Lyapunov exponent measured by other methods 
|nj. By the observation of the 7(^3), its behavior indicate the chaotic transition because is 
directly related with the maximum Lyapunov exponent besides an afin mapping. This is in 
agreement with Eckman et al. [JTJ] statements. Moreover, the algorithms developed and used 
in this work allow us a systematic RQA analysis efficiently. 
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